%% Solve for the steady state
%  and in response to MIT shock to TFP

% asymmetric labor adjustment cost
% fixed mass of potential entrants, each of which gets a signal about 
% their future productivity

% entrants have lower average productivity than do incumbents
% and they slowly grow

% The NICK BLOOM calibration

%- Standard deviation of sales & employment growth
%- autocorrelation of each


%%
addpath ../../../CEtools64/


load eq_mit_mass.mat
eq_mit_kimball = eq_mit;
clear eq_mit;


mit_flag = 3;

%% Set Options
% All
options.Nbell       = 100;        % Number of Bellman iterations before Newton.
options.Nnewt       = 15;       % Maximum number of Newton steps
options.tolc        = 1e-12;     % Tol on value functions
options.T_irf       = 50;       % Number of periods for IRFs
options.T           = 300;     % Number of periods for MIT shock
options.Terg        = 500;      % Number of periods to burn
options.hpisteps    = 100;

% Stationary
options.L = [];
options.itermaxL    = 5000;     % Max number of iterations to find L
options.tolL        = 1e-12;    % Tol for L
options.tolK        = 1e-5;     % Tol for equilibrium K
options.itermaxK    = 100;      % Max iterations for bisection
options.cresult = [];

options.tolD        = 1e-08;
options.itermaxp    = 25;

options.damp = 0;
options.damp_mit = .8;
options.damp_DC  = .9;

optset('bisect', 'tol', 1e-12)

%% Set globals and parameters
% % try vavra numbers
glob.n          = [100,30];   % Number of nodes for b,a
glob.nf         = [100,30];

glob.spliorder  = [1,1];    % Order of splines (Envelope Condition Method seems only robust when quadratic in k)
glob.Na         = 30;           % Number of nodes for quadrature
glob.curv       = .1;            % Amount of curvature for p grid

glob.Nsignal    = 1e6;

glob.minb = 1e-5;
glob.maxb = 50;


%% Model parameters

% preferences
param.beta      = 0.96;               

%  rho_s    phi_L  sigma epsilon   ce_mu     d_E
% 0.60338 0.076289 2.9199   2.835 -2.3445 0.59226

% idiosyncratic shock - set 'em and forget 'em

% note - berger just chooses rho - .81 and sigma = .38 (annual)

%%%%%%%%%%%%% CALIBRATE HERE %%%%%%%%%%%%%
param.sigma_s         = 0.15; % concentration & variance
param.rho_s           = 0.85;  % Concentration & variance
 
% labor adjustment cost
param.phi_L           = 0.05; % as a start, no adjustment costs
 
% Demand parameters
param.sigma           = 50;              % FIX THIS NUMBER, let go of average markup
param.epsilon         = .6731*param.sigma;     % epsilon/sigma is superelasticity, 0 is CES 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

param.omega           = 1;                  
         
param.A = 1;

% mass of operating firms
param.M        =     1.0283; % the mass in the free entry model


param.nu      = .5; % inverse Frisch elasticity
                     % same as in clementi
                     % can calibrate aggregate shock size to hit moments

% free labor adjustment 
param.delta    = 0.19;

glob.mink       = 0;
glob.maxk       = param.sigma^(param.sigma/param.epsilon);

%% Find initial steady state
options.itermax_DC = 100000;

[param,glob]    = setup(param,glob,options);

eq_SS           = solve_eq_ss(param, glob, options);
save('eq_SS.mat', 'eq_SS', 'param', 'glob', 'options');
% param.mu_entry  = eq_SS.mu_entry;
% param.mu_f      = param.mu_entry - param.mu_f_diff*abs(param.mu_entry);

glob.agg = kimball_agg(1, eq_SS, param, glob); % Fix this value!
L        = sum(eq_SS.l.*eq_SS.L);
param.psi = 1/(L^param.nu); % psi = 1/(L^nu*C)

% check that C = 1... 
tosolve = @(C) kimball_agg(C, eq_SS, param, glob) - glob.agg;

fprintf('C is %.5f \n', bisect(tosolve, .5, 1.5));

% %% nvm
% calib = [0.093017 10.753 1.7897 -5.0936  1.1223 2.196 0.62235];
% mom = compute_moments_for_calibration(calib, param, glob, options);

%%
rng(47)
mom = compute_moments(eq_SS, param, glob, options);

%%
disp('------------------------------------')
disp('Targeted Moments')
disp('------------------------------------')
disp(sprintf('Labor regression, within: %f, Target: 0.55', mom.beta_l))
disp('------------------------------------')
disp(sprintf('Hiring autocorr: %f, Target: 0.1281', mom.demp_auto))
disp(sprintf('Sales growth autocorr: %f, target: 0.122', mom.dsales_auto))
disp(sprintf('Emp growth var: %f, Target: 0.0617', mom.labor_churn))
disp(sprintf('Sales growth var: %f, Target: 0.1415', mom.sales_churn))
disp('------------------------------------')
disp(sprintf('Cost weighted markup: %f, Target: 1.25', mom.mu_cw))
disp(sprintf('-------- CONCENTRATION MEASURES -------'))
disp(sprintf('Top 10 percent share: %f, Tagret: 0.75', mom.top10))
disp(sprintf('Frac rel sales below 1: %f, Tagret: 0.75', mom.frac_rel_sales_below_1))
disp('------------------')
%%
addmom = 1;
if addmom == 1
    disp('------------------')
    disp('Additional moments')
    disp('------------------')
    disp(sprintf('JC/JD > 10, big firms: %f', mom.large_jcjd10))
    disp(sprintf('Labor churn: %f', mom.labor_churn))
    disp(sprintf('Sales/labor churn ratio: %f', mom.sales_churn/mom.labor_churn))
    disp(sprintf('Cost size: %f', mom.cost_size))
end

%% Entry shock


if mit_flag == 1
    sim    = set_up_mit_entry(eq_SS, param, glob, options, eq_mit_kimball);
    eq_mit = solve_eq_mit(sim, eq_SS, param, glob, options);
    save('eq_mit_entry_cost.mat'); 
end
%% Exit shock
if mit_flag == 2
    sim    = set_up_mit_exit(eq_SS, param, glob, options);
    eq_mit = solve_eq_mit(sim, eq_SS, param, glob, options);
    save('eq_mit_exit.mat'); 
end
 %% TFP shock
if mit_flag == 3
    sim    = set_up_mit(eq_SS, param, glob, options);
    eq_mit = solve_eq_mit(sim, eq_SS, param, glob, options);
    system('rm eq_mit_tfp.mat'); 
    save('eq_mit_tfp.mat'); 
end
 %% Entry mass shock
if mit_flag == 4
    sim    = set_up_mit_entry_mass(eq_SS, param, glob, options, eq_mit_kimball);
    eq_mit = solve_eq_mit(sim, eq_SS, param, glob, options);
    
    clear eq_mit_ces
    
    save('eq_mit_mass.mat'); 
end

%%
if mit_flag < 1 || mit_flag > 4
    disp('no mit shock')
end